

setwd("C:/Users/Administrator/Desktop/MPhil of Criminology Cambridge/SSRMC lectures/Replication Workshop/Article")

data <- read.csv("Laws of War & Public Opinion_Data_16Jan15.csv", sep=",", header=TRUE)

#check mean of control 1
mean(subset(data,treat1_dummy2=="1. Control",select = dv_complete)$dv_complete,na.rm=TRUE)

#### Creating DV Complete  #####
##transform categorical dependent variables (DV) into 7 point scale

data$dv_complete[data$treatagree==1] <- 7
data$dv_complete[data$treatagree==2] <- 6
data$dv_complete[data$treatneith==1] <- 5
data$dv_complete[data$treatneith==3] <- 4
data$dv_complete[data$treatneith==2] <- 3
data$dv_complete[data$treatdis==2] <- 2
data$dv_complete[data$treatdis==1] <- 1


#create Binary DV (Agree/Disagree)

data$dv_dummy[data$dv_complete>=5] <-1
data$dv_dummy[data$dv_complete<=3] <-0


##Transform data into groups experiencing 8 treatment (IV) conditions

data$group1 <-c(data$treat1_dummy==1 & data$treat2_dummy==1)
data$group2 <-c(data$treat1_dummy==2 & data$treat2_dummy==1)
data$group3 <-c(data$treat1_dummy==3 & data$treat2_dummy==1)
data$group4 <-c(data$treat1_dummy==4 & data$treat2_dummy==1)
data$group5 <-c(data$treat1_dummy==1 & data$treat2_dummy==2)
data$group6 <-c(data$treat1_dummy==2 & data$treat2_dummy==2)
data$group7 <-c(data$treat1_dummy==3 & data$treat2_dummy==2)
data$group8 <-c(data$treat1_dummy==4 & data$treat2_dummy==2)

###Transform Party Variables###
##1.Transform to scale (not used in analysis)##

data$party_complete[data$dem==1] <- 7
data$party_complete[data$dem==2] <- 6
data$party_complete[data$indp==5] <-5 
data$party_complete[data$indp==6] <- 4
data$party_complete[data$indp==4] <- 3
data$party_complete[data$repub==2] <- 2
data$party_complete[data$repub==1] <- 1

##2.Transform to Binary party variables##

data_democrats <- data[which(data$party_complete>5),]
data_republicans <- data[which(data$party_complete<3),]

data_demo1 <- data[which(data$party_complete>=5),]
data_rep1 <- data[which(data$party_complete<=3),]

###Transform demographic variables###

data$male[data$gender==1] <- 1
data$male[data$gender!=1] <- 0

data$race[data$race_1==1] <- "White"
data$race[data$race_2==1] <- "Black"
data$race[data$race_3==1] <- "Asian"
data$race[data$race_4==1] <- "Hispanic"
data$race[data$race_5==1] <- "Other"
data$race[data$race_6==1] <- "Other"
data$race[data$race_7==1] <- "Other"

data$College[data$educ_1==1] <- 0
data$College[data$educ_2==1] <- 0
data$College[data$educ_3==1] <- 0
data$College[data$educ_4==1] <- 0
data$College[data$educ_5==1] <- 1
data$College[data$educ_6==1] <- 1
data$College[data$educ_7==1] <- 1
data$College[data$educ_8==1] <- 1
data$College[data$educ_9==1] <- 1
data$College[data$educ_10==1] <- 1


#################
#    Results    #
#################

####### Figure 1: Summary Statistics #######

summary(data$male) ## why?

summary(data$Age,na.rm=TRUE)
sd(data$Age,na.rm=TRUE)

summary(data$College,na.rm=TRUE)
sd(data$College,na.rm=TRUE)

## republicans ##

data$rep[data$repub==1 & data$repub==2] <- 1
data$rep[data$indp==4] <- 0
data$rep[data$indp==5] <- 0
data$rep[data$indp==6] <- 0
data$rep[data$dem==1] <- 0
data$rep[data$dem==2] <- 0
summary(data$rep)
sd(data$rep)

data$citizen2[data$citizen==1] <- 1
data$citizen2[data$citizen!=1] <- 0
summary(data$citizen2)

data$White[data$race=="White"] <- 1
data$White[data$race!="White"] <- 0
summary(data$White,na.rm=TRUE)

data$Black[data$race=="Black"] <- 1
data$Black[data$race!="Black"] <- 0
summary(data$Black,na.rm=TRUE)

## the same to Asian, Hispanic, and Other ##


  ##########################################################
 #    Table 3: Mean Response and Confidence Intervals     #
##########################################################

### control 1
mean(subset(data,group1,select = dv_complete)$dv_complete,na.rm=TRUE) ## or
t.test(subset(data,group1,select = dv_complete)$dv_complete,na.rm=TRUE,conf.level =.95)

### control 2
mean(subset(data,group2,select = dv_complete)$dv_complete,na.rm=TRUE) ## or
t.test(subset(data,group2,select = dv_complete)$dv_complete,na.rm=TRUE,conf.level =.95)

### control 3
mean(subset(data,group3,select = dv_complete)$dv_complete,na.rm=TRUE) ## or
t.test(subset(data,group3,select = dv_complete)$dv_complete,na.rm=TRUE,conf.level =.95)

### control 4
mean(subset(data,group4,select = dv_complete)$dv_complete,na.rm=TRUE) ## or
t.test(subset(data,group4,select = dv_complete)$dv_complete,na.rm=TRUE,conf.level =.95)

### control 5,6,7,8
t.test(subset(data,group5,select = dv_complete)$dv_complete,na.rm=TRUE,conf.level =.95)
t.test(subset(data,group6,select = dv_complete)$dv_complete,na.rm=TRUE,conf.level =.95)
t.test(subset(data,group7,select = dv_complete)$dv_complete,na.rm=TRUE,conf.level =.95)
t.test(subset(data,group8,select = dv_complete)$dv_complete,na.rm=TRUE,conf.level =.95)

### total1234 & 5678
t.test(subset(data,treat2_dummy==1, select=dv_complete)$dv_complete,na.rm=TRUE,conf.level =.95)
t.test(subset(data,treat2_dummy==2, select=dv_complete)$dv_complete,na.rm=TRUE,conf.level =.95)

### total15, 26, 37, 48
t.test(subset(data,treat1_dummy==1, select=dv_complete)$dv_complete,na.rm=TRUE,conf.level =.95)
t.test(subset(data,treat1_dummy==2, select=dv_complete)$dv_complete,na.rm=TRUE,conf.level =.95)
t.test(subset(data,treat1_dummy==3, select=dv_complete)$dv_complete,na.rm=TRUE,conf.level =.95)
t.test(subset(data,treat1_dummy==4, select=dv_complete)$dv_complete,na.rm=TRUE,conf.level =.95)

### total all
t.test(data$dv_complete,na.rm=TRUE,conf.level =.95)


  ##############################################################
 #####     Figure 1: The Effect of International Law     ######
##############################################################

install.packages("sciplot")
library("sciplot")

data_1 <- data[which(data$treat2_dummy==1),]
data_2 <- data_1[which(data_1$treat1_dummy<3),]

lineplot.CI(x.factor = treat1_dummy2, response = dv_complete, data =  data_2,cex = 2,ylab = "Approve of Continued Bombing", xlab="Treatment Group",cex.lab = 1,x.leg = 1,col = c("black","black"), pch = c(16,16), ylim=c(2,4), lty=0, fun = function(x) mean(x, na.rm=T), ci.fun = function(x) c(mean(x, na.rm=T)-1.96*se(x, na.rm=T), mean(x, na.rm=T)+1.96*se(x, na.rm=T)))

mean(subset(data,group1,select = dv_dummy)$dv_dummy,na.rm=TRUE)- mean(subset(data,group2,select = dv_dummy)$dv_dummy,na.rm=TRUE)

t.test(subset(data,group1,select = dv_complete)$dv_complete,na.rm=TRUE,subset(data,group2,select = dv_complete)$dv_complete,na.rm=TRUE,conf.level =.95)

  ##################################################################
 ######  Figure 2: The Effect of IL & Morality Treatments   #######
##################################################################

### subsetting data
data_1 <- data[which(data$treat2_dummy==1),]

lineplot.CI(x.factor = data$treat1_dummy2, response = data$dv_complete, data =  data_1,cex = 2, ylab = "Approve of Continued Bombing", xlab="Treatment Group",cex.lab = 1, x.leg = 1, col = c("black","black"), pch = c(16,16), ylim=c(2,4),lty=0, fun = function(x) mean(x, na.rm=T), ci.fun = function(x) c(mean(x, na.rm=T)-1.96*se(x, na.rm=T), mean(x, na.rm=T)+1.96*se(x, na.rm=T)))

mean(subset(data,group1,select = dv_dummy)$dv_dummy,na.rm=TRUE)- mean(subset(data,group3,select = dv_dummy)$dv_dummy,na.rm=TRUE)
t.test(subset(data,group1,select = dv_complete)$dv_complete,na.rm=TRUE,subset(data,group3,select = dv_complete)$dv_complete,na.rm=TRUE,conf.level =.95)

mean(subset(data,group3,select = dv_dummy)$dv_dummy,na.rm=TRUE)- mean(subset(data,group2,select = dv_dummy)$dv_dummy,na.rm=TRUE)
t.test(subset(data,group2,select = dv_complete)$dv_complete,na.rm=TRUE,subset(data,group3,select = dv_complete)$dv_complete,na.rm=TRUE,conf.level =.95)

mean(subset(data,group4,select = dv_dummy)$dv_dummy,na.rm=TRUE)- mean(subset(data,group2,select = dv_dummy)$dv_dummy,na.rm=TRUE)
t.test(subset(data,group2,select = dv_complete)$dv_complete,na.rm=TRUE,subset(data,group4,select = dv_complete)$dv_complete,na.rm=TRUE,conf.level =.95)

mean(subset(data,group3,select = dv_dummy)$dv_dummy,na.rm=TRUE)- mean(subset(data,group4,select = dv_dummy)$dv_dummy,na.rm=TRUE)
t.test(subset(data,group3,select = dv_complete)$dv_complete,na.rm=TRUE,subset(data,group4,select = dv_complete)$dv_complete,na.rm=TRUE,conf.level =.95)



  #############################################################
 ######     Figure 3: The Effect of Reciprocity     ##########
#############################################################

mean(data_1$treat1_dummy2)
mean(data_1$dv_complete,na.rm=T)

lineplot.CI(x.factor = data$treat1_dummy2, response = data$dv_complete, data =  data_1,cex = 2, ylab = "Approve of Continued Bombing", xlab="No Reciprocity",cex.lab = 1, x.leg = 1, col = c("black","black"), pch = c(16,16), ylim=c(2,4),lty=0, fun = function(x) mean(x, na.rm=T), ci.fun = function(x) c(mean(x, na.rm=T)-1.96*se(x, na.rm=T), mean(x, na.rm=T)+1.96*se(x, na.rm=T)))

##### subset group with reciprocity treatment

data_3 <-data[which(data$treat2_dummy==2),]
            
lineplot.CI(x.factor = data$treat1_dummy3, response = data$dv_complete, data =  data_3,cex = 2, ylab = "Approve of Continued Bombing", xlab="Reciprocity",cex.lab = 1, x.leg = 1, col = c("black","black"), pch = c(16,16), ylim=c(2,4),lty=0, fun = function(x) mean(x, na.rm=T), ci.fun = function(x) c(mean(x, na.rm=T)-1.96*se(x, na.rm=T), mean(x, na.rm=T)+1.96*se(x, na.rm=T)))

mean(data_1$dv_dummy,na.rm=T)- mean(data_3$dv_dummy,na.rm=TRUE)
wilcox.test(data_1$dv_complete,data_3$dv_complete)

  #####################################################################################
 ######     Figure 4: International Law's Effect by Political Affiliation    #########
#####################################################################################

#### subset data with political affiliation and international law treatment
data_demIL <- data_democrats[which(data_democrats$treat1_dummy<3),]
data_repIL <- data_republicans[which(data_republicans$treat1_dummy<3),]

##### run models

lineplot.CI(x.factor = data_demIL$treat1_dummy2, response = data_demIL$dv_complete, data =  data_demIL,cex = 2, ylab = "Approve of Continued Bombing", xlab="Treatment Group",cex.lab = 1, x.leg = 1, col = c("black","black"), pch = c(16,16), ylim=c(2,5),lty=0, main="Democrats", fun = function(x) mean(x, na.rm=T), ci.fun = function(x) c(mean(x, na.rm=T)-1.96*se(x, na.rm=T), mean(x, na.rm=T)+1.96*se(x, na.rm=T)))

lineplot.CI(x.factor = data_repIL$treat1_dummy2, response = data_repIL$dv_complete, data =  data_repIL,cex = 2, ylab = "Approve of Continued Bombing", xlab="Treatment Group",cex.lab = 1, x.leg = 1, col = c("black","black"), pch = c(16,16), ylim=c(2,5),lty=0, main="Republicans", fun = function(x) mean(x, na.rm=T), ci.fun = function(x) c(mean(x, na.rm=T)-1.96*se(x, na.rm=T), mean(x, na.rm=T)+1.96*se(x, na.rm=T)))

###### difference between republicans and democrats

wilcox.test(subset(data_demIL,treat1_dummy==1,select = dv_complete)$dv_complete,subset(data_demIL,treat1_dummy==2,select = dv_complete)$dv_complete)

wilcox.test(subset(data_repIL,treat1_dummy==1,select = dv_complete)$dv_complete,subset(data_repIL,treat1_dummy==2,select = dv_complete)$dv_complete)

t.test(subset(data_repIL,treat1_dummy==1,select = dv_complete)$dv_complete,subset(data_repIL,treat1_dummy==2,select = dv_complete)$dv_complete)

####### difference: following the original coding scheme for political affiliation

data_demIL1 <- data_demo1[which(data_demo1$treat1_dummy<3),]
data_repIL1 <- data_rep1[which(data_rep1$treat1_dummy<3),]

wilcox.test(subset(data_demIL1,treat1_dummy==1,select = dv_complete)$dv_complete,subset(data_demIL1,treat1_dummy==2,select = dv_complete)$dv_complete)

wilcox.test(subset(data_repIL1,treat1_dummy==1,select = dv_complete)$dv_complete,subset(data_repIL1,treat1_dummy==2,select = dv_complete)$dv_complete)
t.test(subset(data_repIL1,treat1_dummy==1,select = dv_complete)$dv_complete,subset(data_repIL1,treat1_dummy==2,select = dv_complete)$dv_complete)

####################################
#####   Improvements   #######  ### Wilcoxon, effect size, missingness (see email), multi-inputation ###
##################################

############################################
####  Check approriateness for T test  #####
############################################

###### Check Normal Distribution for using the T test ##########

qqnorm(data_demIL$dv_complete)
hist(data_demIL$dv_complete)
hist(data_repIL$dv_complete)
hist(data_1$dv_complete)
hist(subset(data,group1,select = dv_complete)$dv_complete)
hist(subset(data,group2,select = dv_complete)$dv_complete)
hist(subset(data,group3,select = dv_complete)$dv_complete)
hist(subset(data,group4,select = dv_complete)$dv_complete)
hist(subset(data,group5,select = dv_complete)$dv_complete)
hist(subset(data,group6,select = dv_complete)$dv_complete)
hist(subset(data,group7,select = dv_complete)$dv_complete)
hist(subset(data,group8,select = dv_complete)$dv_complete)
hist(data_2$dv_complete)
hist(data_3$dv_complete)


###### Apply Mann�CWhitney test rather than T-test: because data is not normally distributed #######

### group 1, group2
wilcox.test(subset(data,group1,select = dv_complete)$dv_complete,subset(data,group2,select = dv_complete)$dv_complete)
t.test(subset(data,group1,select = dv_complete)$dv_complete,subset(data,group2,select = dv_complete)$dv_complete,na.rm=T)

### group2, group3
wilcox.test(subset(data,group2,select = dv_complete)$dv_complete,subset(data,group3,select = dv_complete)$dv_complete)
t.test(subset(data,group2,select = dv_complete)$dv_complete,subset(data,group3,select = dv_complete)$dv_complete,na.rm=T)

wilcox.test(subset(data,group1,select = dv_complete)$dv_complete,subset(data,group3,select = dv_complete)$dv_complete)
wilcox.test(subset(data,group2,select = dv_complete)$dv_complete,subset(data,group4,select = dv_complete)$dv_complete)
wilcox.test(subset(data,group3,select = dv_complete)$dv_complete,subset(data,group4,select = dv_complete)$dv_complete)

wilcox.test(data_1$dv_complete,data_3$dv_complete)
mean(data_1$dv_dummy,na.rm=T) - mean (data_3$dv_dummy,na.rm=T)

wilcox.test(subset(data,group1,select = dv_complete)$dv_complete,subset(data,group5,select = dv_complete)$dv_complete)
wilcox.test(subset(data,group2,select = dv_complete)$dv_complete,subset(data,group5,select = dv_complete)$dv_complete)
wilcox.test(subset(data,group3,select = dv_complete)$dv_complete,subset(data,group5,select = dv_complete)$dv_complete)
wilcox.test(subset(data,group4,select = dv_complete)$dv_complete,subset(data,group5,select = dv_complete)$dv_complete)

wilcox.test(subset(data,group1,select = dv_complete)$dv_complete,subset(data,group6,select = dv_complete)$dv_complete)
wilcox.test(subset(data,group2,select = dv_complete)$dv_complete,subset(data,group6,select = dv_complete)$dv_complete)
wilcox.test(subset(data,group3,select = dv_complete)$dv_complete,subset(data,group6,select = dv_complete)$dv_complete)
wilcox.test(subset(data,group4,select = dv_complete)$dv_complete,subset(data,group6,select = dv_complete)$dv_complete)

wilcox.test(subset(data,group1,select = dv_complete)$dv_complete,subset(data,group7,select = dv_complete)$dv_complete)
wilcox.test(subset(data,group2,select = dv_complete)$dv_complete,subset(data,group7,select = dv_complete)$dv_complete)
wilcox.test(subset(data,group3,select = dv_complete)$dv_complete,subset(data,group7,select = dv_complete)$dv_complete)
wilcox.test(subset(data,group4,select = dv_complete)$dv_complete,subset(data,group7,select = dv_complete)$dv_complete)

wilcox.test(subset(data,group1,select = dv_complete)$dv_complete,subset(data,group8,select = dv_complete)$dv_complete)
wilcox.test(subset(data,group4,select = dv_complete)$dv_complete,subset(data,group8,select = dv_complete)$dv_complete)
wilcox.test(subset(data,group3,select = dv_complete)$dv_complete,subset(data,group8,select = dv_complete)$dv_complete)
wilcox.test(subset(data,group2,select = dv_complete)$dv_complete,subset(data,group8,select = dv_complete)$dv_complete)
mean(subset(data,group3,select = dv_dummy)$dv_dummy,na.rm=T) - mean(subset(data,group8,select = dv_dummy)$dv_dummy,na.rm=T)

wilcox.test(subset(data,group5,select = dv_complete)$dv_complete,subset(data,group1,select = dv_complete)$dv_complete)
mean(subset(data,group1,select = dv_dummy)$dv_dummy,na.rm=T) - mean(subset(data,group5,select = dv_dummy)$dv_dummy,na.rm=T)


###### Check the number of missing value #########

summary(data$treatagree)
is.na(data$treatagree)

sum(is.na(data$treatagree))/length(data$treatagree)

